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ABSTRACT 

Interfaces between hot and cold magnetized plasmas exist in various astro- 
physical contexts, for example where hot outflows impinge on an ambient in- 
terstellar medium (ISM). It is of interest to understand how the structure of 
the magnetic field spanning the interface affects the temporal evolution of the 
temperature gradient. Here we explore the relation between the magnetic field 
topology and the heat transfer rate by adding various fractions of tangled vs. 
ordered field across a hot-cold interface allow the system to evolve to a steady 
state. We find a simple mathematical relation for the rate of heat conduction as 
a function of the initial ratio of ordered to tangled field across the interface. We 
discuss potential implications for the astrophysical context of magnetized wind 
blown bubbles (WBB) around evolved stars. 

Subject headings: magneto- hydrodynamics, planetary nebula, magnetic reconnection, 
wind blown bubbles, anisotropic heat conduction 
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Introduction 



Interfaces between hot and cold plasmas can occur in astrophysics where understanding 
the rate of thermal conduction may be an important part of the the astrophysical 
phenomenology. One example occurs in wind blown bubbles (WBB) of evolved stars where 
magnetized hot supersonic outflow shock heats the cooler ambient magnetized interstellar 
medium. For such WBB, there are examples where the presumed shock heated bubble 



is cooler than expected if only radiative cooling is considered (Zhekov et al. (2011)). A 



possible explanation is that heat loss through the interface of hot bubble into the cold 



shell via thermal conduction reduces the temperature of the hot bubble (Zhekov & Park 



(1998), Zhekov & Myasnikov (2000)). However the source of heat into the cold side of 
the interface will continuously evaporate material there and potentially induce interface 
instabilities and mass mixing (Stone & Zweibel ( 2009[ )) that could tangle the magnetic 
field. Understanding the thermal conduction and its dependence on magnetic structure 
is important for determining the thermal properties of the plasma on either side of the 
interface. 

A second example is the unexpected slow mass deposition rate of the cooling flows in 



some galaxy cores which might be inhibited by a restricted thermal conduction (Rosner 



fc Tucker] ( |1989D , [Balbus fc Reynolds] ( |2008D , [MikeUides et al.| ( |201lD ). In the intracluster 
medium (ICM), the tangled magnetic field can potentially produce a strongly anistropic 
thermal conductivity that may significantly influence temperature and density profiles 



Chandran & Maron 


(2004 


); 


Maron et al. 


(2004 


); 



et al. (2011)) 



For the ISM and ICM, it is usually valid to assume that the electrons are totally 
inhibited from moving across field lines (McCourt et al ( 2011[ )), as the electron mean free 
path is much greater than the electron gyroradius. The magnetic field structure therefore 
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plays a key role in controlling the rate of thermal conduction since electrons can move freely 
only along the field lines. The result is a strong thermal conductivity parallel to the field 
lines and a weak conductivity across the field lines. 

The quantitative subtleties of how a complicated magnetic field structure affects 
thermal conduction for raises the open question of whether there is a simple measure of 
field tangling that allows a practical but reasonably accurate correction to the isotropic 
conduction coefficient for arbitrarily tangled fields. In this context, two classes of problems 
can be distinguished. The first is the conduction in a medium for which forced velocity flows 
drive turbulence, which in turn tangles the field into a statistically steady state turbulent 



spectrum (Tribble (1989), Tao (1995); Maron et al. (2004)). The second is the case in 



which the flow is laminar and the level of conduction inhibition is compared when the field 
starts from initial states of different levels of tangling subject to an imposed temperature 
difference across an interface. This second problem is the focus of our preset paper. 

Using the ASTROBEAR magnetohydrodynamics code with anisotropic thermal 
conduction, we investigate the influence of initial magnetic structure on thermal conduction 
in an otherwise laminar flow. The key questions we address are: (1) does the interface 
become unstable? (2) how fast is the thermal conduction across the interface compared to 
the unmagnetized case? 

We study these questions using different initial magnetic configurations imposed on a 
planar thermal interface to determine how the conduction depends on the amount of field 
tangling across the interface. 

In section 2, we review the basic equations of MHD with anisotropic thermal 
conduction. In sections 3 and 4 we provide detailed description of the simualtion setup. In 
section 5 and 6 we present the simulation results and analyses. In section 7, we discuss the 
simulation results in the context of the WBB cooling problem and the cooling flow problem 
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in cores of galaxy clusters. The appendix provides more detailed information on the testing 
of the ASTROBEAR code. 



2. MHD Equations with Anisotropic Heat Conduction 

The MHD equations with anisotropic heat conduction that we will solve are given by: 

|^ + V-(pv) = 0, (1) 

g^ + Vx(vxB) = 0, (3) 

f +V.[v(B + p+f!)-5^f:^l + V.Q = 0, (4) 

where p, v, B and p are the density, velocity, magnetic field, and pressure, and E denotes 
the total energy given by 

vv B B 



where the internal energy e is given by 



and 7 = 5/3. In our simulations, we will assume that the heat flux is confined to be 
parallel to the magnetic field lines. This assumption applies only when the ratio of electron 
gyro-radius to field gradient scales is small. Under this assumption, the heat flux parallel 
to field lines can be written as 

g = -/.ii(vr)||, (7) 

where the subscript || indicates parallel to the magnetic field, and k\\ is the classical Spitzer 
heat conductivity: k\\ — KcT"^'^, with aCc = 2 x cm s K~'^-^ . We take k\\ to be a 
constant throughout our simulations and so hereafter write it simply as k. 
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The ASTROBEAR code uses an operator splitting method to solve these MHD 
equations with heat conduction. The viscosity and resistivity are ignored in our calculation 
so the dissipation is numerical only. The ideal MHD equations are then solved with the 
MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) primitive method 
with TVD (Total Variation Diminishing) preserving Runge-Kutta temporal interpolation. 
The result is then sent to the implicit linear solver utilizing the High Performance 
Preconditioners (Hypre) to solve the anisotropic heat conduction equation. The linear 
solver requires temporal sub-cycling technique to maintain its accuracy. The code runs in 
parallel with fixed grid domain. 



3. Problem Description and Analytical Model 

Our initial set up involves hot and cold regions separated by a thin planar interface. 
We study how the magnetic field configuration alters the heat transfer rate between the hot 
and cold regions in presence of anisotropic heat conduction. We study the problem in 2-D. 

To guide subsequent interpretation of the results, we first compare two simple but 
illustrative limits of magnetic field orientation: (1) a uniform magnetic field aligned with 
the direction normal to the interface; (2) a uniform magnetic field perpendicular to the 
normal direction of the interface. In case (1), because the angle between the magnetic 
field and temperature gradient is everywhere zero, heat conduction across the interface is 
expected to take on the Spitzer value associated with isotropic heat conduction. In case (2) 
however, the angle between the magnetic field and the temperature gradient is always 90°, 
so with our approximations, heat cannot fiow across the interface. 

We define a heat transfer efficiency ( equal to the magnetic field-regulated heat transfer 
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rate divided by the isotropic Spitzer rate, namely, 




where q is defined as the amount of thermal energy transported through the interface per 
unit time. The average angle 9 between the temperature gradient and the uniform magnetic 
field then plays an important role in determining C,. At = 0, C = 1. At = 7r/2, C = 0- 

We now address the infiuence of both a mean field and a tangled field Consider 
there to be a strongly tangled local field that has no mean value in the direction normal to 
the interface, i.e. Bo,a; whose total magnitude is Bq, and a global magnetic field aligned 
with the normal of the interface of magnitude 3^. If 3^ ^ i?o, the magnetic field around 
the interface only slightly deviates from the normal direction and C should be close to 1. If 
Bd <^ Bq, one would expect that the global energy transfer would be slow and C should be 
close to zero. 

If Bd and Bq are comparable, we expect < C < 1. We also expect C can change 
throughout the evolution if the strucuture of the magnetic field is modified by the dynamics 
of heat transfer. It is instructive to ask whether the feedback from the magnetic field 
structure evolution will amplify the heat transfer by creating more channels, or shut it 
down. The answer depends on the infiuence of magnetic reconnection, as we will see 
from the simulations. Only if magnetic reconnection acts to smooth out local small scale 
structures and link the initially isolated structures to the global mean field across the 
interface then we would expect the heat conductivity to increase. 

EDITOR: PLACE FIGURE □ HERE. 

In what follows, we refer to the initial tangled field region as "the interaction region". 
Figure llVa) shows a schematic of initial and hypothetical evolved steady state field 
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configurations for sucli a tangled field set up. From tlie figure we can see tliat the initial 
field configuration forms a "wall" which restricts energy transfer across the two interaction 
region. However, if the subsequent evolution evolves to the steady state shown in (b), then 
expansion of the interaction region and magnetic reconnection has allowed the field to 
penetrate through the entire region. Thus the initial "wall" of tangle field wall is destroyed 
and thermal conduction will be less inhibited than initially. We will check how accurately 
this proposed picture of destruction of field wall is valid from analyzing our numerical 
simulations, and quantitatively discuss the effects on the energy transfer. 



For our initial conditions, we set up an interface between hot and cold regions in 
mutual pressure equilibrium. The temperature distribution on the horizontal (x) axis is 
given by: 



in the region < x < 1 with Tq = 100 in computational units. This temperature profile 
has a sharp gradient at x = 0. The temperature distribution is plotted in Figure [2]^a). The 
region 0.4 < x < 0.5 is the interaction region we described in the previous section. At the 
two side boundaries, the temperature is set to be constant, and uniform across the regions 
of each respective side of the box connecting to that side of the interaction region. We are 
primarily interested in the region of the box where the heat transfer occurs and noticeably 
evolves during the simulation run time. This means we will mainly focus on the interaction 
region. The horizontal length of the interaction region in the simulation domain is 0.1 in 
computational units. 



4. Simulation Setup 




(9) 
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The thermal pressure is set to be in equihbrium over the entire box, that is 

P{x) = Po (10) 
with Po = 100. The density distribution is set up by the ideal gas law, namely: 



in computational units. 



For the Spitzer diffusion coefficient, we assume the diffusion is linear as in Eq.(7), and 
use the approximation: /ty = HcT^f^, where Kc is the classical conductivity, and Tmid is 
taken to be the middle value of temperature across the interface, about 0.5 Tq. 

We choose the initial field configuration: 

B, = Ba + Bosm{mry/X), (12) 

By = Bq sm{mi X / \) (13) 

where n and A are the mode number and wavelength of the tangled field respectively, 
Bq = 10~^ in computational units, and 5^ can assume various initial values that reflect the 
evolving global field as the result of reconnection. This initial field configuration is therefore 
one of a locally tangled field surrounding the interface with one measure of the tangle given 
by: 

R = Bd/Bo (14) 

When R = 0, there are only locally confined field lines, whereas R = oo indicates a 
straight horizontal field without any "tangling". As R increases, the relative fraction of 
field energy corresponding to lines which penetrate through the interaction region increases. 
In our simulations, we consider cases with R = 0.0, 0.2, 0.4, 0.6, 1, 2, 4, oo. Figure |2](a). 
Figure |4]^a) and Figure [5|a) show the magnetic field configuration for initial R values of 0.0, 
0.4, 1.0. 
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We note that our MHD approximation a priori implies that the electron gydroradius is 
much smaller than the length scale of one grid cell. Thus the dissipation scale and all field 
gradient scales are larger than the electron gyro-radius by construction in our simulations. 

We run simulations with typical resolution of 2048 cells on the horizontal axis in fixed 
grid mode. Runs with doubled resolution showed no significant differences compared to 
the standard resolution runs. We use fixed boundary conditions at the x boundaries: the 
pressure, density and temperature at the two ends are fixed to their initial values, as is the 
magnetic field. We use periodic boundary conditions for the y-axis boundaries. 

There are five parameters whose infiuence determine the simulation behavior and guide 
interpretation of results: 

1. Plasma /3. /3 = ^0- has little effect on diffusion because even with very high values 
of the plasma j3 used in the simulation, we are still in the MHD regime and the gyro-radii 
of electrons are assumed small. Thus the direction of thermal conduction is not locally 
affected by /3. It is possible that instabilities could arise in the low /3 limit that affect 
pressure balance during the evolution of the simulations but that turns out not to be the 
case for the /3 range of 10^ ~ 10* that we use. The value of /3 in this range does not exihbit 
any infiuence on the simulation result as indicated by our numerical experiment. 

2. Initial Tangle measure R = 3^/ Bq. If R » 1, the local small scale field can 
mostly be ignored and Spitzer thermal conductivity is expected, whereas if i? << 1 a value 
much less than Spitzer is expected. 

3. Ratio of the diffusion time scale to the sound crossing time scale for one 
grid cell: 



(15) 
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where p is the density, I is the characteristic gradient length scale of temperature: 
I = min{j^-jr^) and Cg is the sound speed. If r << 1, thermal diffusion would initially 
dominate and the pressure equilibrium would be broken by this fast energy transfer. If 
r » 1, then the pressure equilibrium would be well maintained throughout the entire 
evolution and the energy transfer may be viewed as a slow relaxation process. In our 
simulation, r ~ 0.3 initially, so that diffusion induces a pressure imbalance. Eventually, as 
the heat transport slows, the pressure equilibrium catches up and is maintained. 

4. Ratio between the temperature gradient scale length and the wavelength 
of the tangled field: h — 2ttI/X — kl . li h — there is no tangled field, and no inhibition 
to heat transfer. As h inreases, the field becomes more tangled, and the energy is harder 
to transfer. However, a large h value may also result in increased magnetic reconnection, 
because the Lundquist number of field confined in a smaller region is larger, for the same 
field strength. Thus would then lower h. 

5. Mean global energy transfer rate: q = SE/tbai, where tbai is defined as the 
time needed for the hot region and cold region to reach a certain degree of temperature 
equilibrium by a transfer of heat energy 5E across the interface. 

A mathematical expression for the heat transfer rate can be derived by considering a 
slab with a planar interface aligned with the y direction at the middle of the interaction 
region (x = 0) with a tangled magnetic field, and an average temperature gradient aligned 
in the x direction. Define the global temperature gradient as |VT|p = {Thot — Tcoid) / {Tq L) , 
where the subscripts "hot" and "cold" denote the characteristic temperatures of the hot 
and cold regions, L is the width of the interaction region, and Tq is a normalization factor 
which ensures |VT|g has dimensions of inverse length. We can then integrate over the 
volume of the interaction region (and since there is no z-dependence, the essential content 
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is an area integral) to obtain the effective fieat flux through this region: 

q^D\WT\g j ^^dxdy, (16) 

where D is a constant that depends on neither the magnetic field nor the temperature 
distribution, |-B| is the local field strength. Notice that this expression is valid only when 
the magnetic field is varying at a length scale smaller than the interaction region length. 

Using Eqs.(12), (13) and (14) in (16) and the approximation that the areal average in 
the interaction region (Bq ■ B^) ~ so that ((Bq + B^^)^) ~ {Bq + B^), we obtain 

q^D^Ik£. (17) 

For the unmagnetized isotropic case, or for transfer with a field entirely aligned with the 
temperature gradient, we have instead 

q,^D\VT\,. (18) 

Dividing Eq.(18) by Eq.(19), we obtain an appoximation for the heat transfer efficiency 
over the interaction region: 

, ^ ■ (19) 

If the initial temperature profiles are identical for different field configurations, this formula 
can then be used to estimate the expected energy tranfer rate from situations with various 
field configuration. By normalizing the heat transfer rate to that of the isotropic heat 
conduction case, we obtain the heat transfer efficiency C,. The accuracy of Eq.(19) can 
be tested by plotting the heat transfer efficiency obtained from the simulations against 
measured values of R. 



If magnetic reconnection occurs during the time evolution of the heat transfer process, 
then conduction channels can open up and the energy exchange can be accelerated. We 
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would then expect the actual curve of ( vs R to evolve to be higher than the value Eq.(19) 
predicts in situations with low R values. Meanwhile, for high R, the analytical prediction 
and the real physical outcome should both approach the horizontal line ( = 1, which 
denotes conductive efficiency consistent with the unmagnetized case. We emphasize that 
R as used in this paper is always calculated with the the initial values of the magnetic 
field, not time evolved values, and that Eq.(19) is valid when estimating a cold to hot 
interface with initial tangle measure as the ratio of initial global straight field to initially 
local tangled field. To follow a measure of the tangle that evolves with time, a generalized 
tangle measure should be calculated in a more sophisticated manner and the integral form 
(Eq.(16)) should be applied. 



5. Simulation Results 

We choose initial conditions with values R = 0.0, 0.2, 0.4, 0.6, 1.0, 2.0, 4.0 to run the 
simulations. The simulation run time is taken to be 1.2 (which corresponds to 12, 000 
years in real units for WBB. The initial cuts of temperature and magnetic field lines 
for R = 0.0, 0.4, 1.0 are shown in Figure [2]^a), Figure |4]^a) and Figure [sj^a) respectively. 
Figure |3](a) shows the initial cut of the density distribution in the R = 0.0 run. We also run 
simulations with purely horizontal magnetic field lines, equivalent to the R = oo case, and 
runs with purely vertical field lines. Frames (b) to (d) in Figure [2] to Figure [s] are from the 
late stages of the evolution, and the final frames always display the steady state of the runs. 
A steady state is facilitated by the fact that the boundaries are kept at a fixed temperature 
throughout the simulations. 



EDITOR: PLACE FIGURE [2] HERE. 
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In Figure [6} we plot the mean cuts of the temperature Tc, obtained by averaging 
the temperature along y axis, against the x position for selected evolution times. Since 
the anisotropic heat conduction is initially faster than the pressure equilibration rate, the 
energy distribution around the temperature interface change rapidly until about t = 0.4. 
This energy transfer is mostly confined to the interaction region for the low Rq runs, since 
in these cases only a few field lines can penetrate into the entire interaction region. 

EDITOR: PLACE FIGURE i HERE. 



EDITOR: PLACE FIGURE i HERE. 



EDITOR: PLACE FIGURE [5] HERE. 

During the initial heat exchange phase, the thermal energy and density quickly 
redistribute in the interaction region. As seen in Figure |2]^b), islands at x = 0.48 are formed 
by material bounded by the magnetic field lines, since the field orientation blocks heat 
exchange with the surroundings. Around x = 0.4, there are also cavities formed where the 
thermal energy is inhibitted from flowing. The magnetic field lines, which form complete 
sets of loops in the R = 0.0 case, begin to distort. It can be observed that the field lines 
are more strongly distorted in the low density part of the interaction region than in the 
high density part. This occurs because velocity gradients are driven by the early rapid 
redistribution of heat (pressure) by conduction. 

At time t = 0.4 (see Figure |2|^b)), the field lines surrounding the cavities at x = 0.4 
reconnect, making thermal exchange possible. During the evolution, field lines begin to link 
the interaction region to the hot material on the left. This phenomenon is most apparent 
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in Figure [2]^d), which marks the final state of the thermal energy exchange. We also see 
that there is little difference between Figure gc) and Fi gure|2](d), because at late stage of 
the process, the thermal diffusion gradually slows so that the magnetic field configuration 
approaches a steady state. 

By comparing Figure |6]^c) with Figure |6](d), we see that the mean cuts of temperature 
show little difference for all values of R. The mean cuts of temperature Tc exhibit a jump in 
the region of x = 0.35 ~ 0.5, but are relatively smooth on either side of this region. This 
shows that even though the tangled field "wall" has been broken and allows channels of 
thermal conduction through it, the temperature profiles is not as smooth as in the purely 
straight field case. 

EDITOR: PLACE FIGURE [6] HERE. 

For the cases of i? = 0.4, there are field lines which penetrate the entire interaction 
region from the start. By observing the evolution of the magnetic field lines at about 
X = 0.38, we see that magnetic reconnection is still happening, and causes the field loops 
to merge. The observed behavior resembles the process displayed by Figure [T] When 
-R = 1, there are hardly any temperature islands that bounded by magnetic field loops. The 
evolution of the field lines shows less dramatic reconnection and evolve in what appears as 
more gentle straightening. 



6. Discussion 

EDITOR: PLACE FIGURE [7] HERE. 



We begin our analysis with the evolution of the heat fiux. The average heat fiux per 
computation cell for different values of R is plotted as a function of time in Figure IWa). 
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Note that in the vertical field case (B = Bylj) the heat flux remains zero as fleld entirely 
inhibits electron motion across the interface. For cases with R > 1, the heat flux decreases 
throughout the evolution. Recall that R > 1 implies cases where the "tangled" portion of 
the fleld is relatively weak and heat is quickly transported from one side of the interface 
to the other. Thus the trend we see for R > 1 occurs as the temperature distribution 
approaches its equilibrium value. For lower R values, especially those of R < 0.5, an initial 
phase of heat flux ampliflcation is observed as magnetic reconnection in the early evolution 
opens up channels for heat to transfer from hot to cold regions. At the late stage of the 
evolution when reconnection has established pathways from deeper within the hot region to 
deeper within the cold region temperature equilibration dominates leading to a decreasing 
heat flux phase as observed in the R> 1 cases. Note that the similarity between the R> 2 
cases and the R = oo case is predicted by Eq.(19): as the global fleld comes to dominate, 
the heat flux inhibition imposed by anisotropic heat conduction in the local tangled fleld 
can be ignored. 

EDITOR: PLACE FIGURE i HERE. 

In order to understand the influence of magnetic reconnection on heat transfer rates 
we compare simulations with different fllling fractions of the tangled fleld. Two cases are 
shown in Figure [sj^a): (1) a temperature interface with a "volume fllling" tangled fleld 
and (2) temperature interface with the tangled fleld fllling only the region surrounding the 
interface. In case (2) the rest of the domain is fllled with straight fleld lines connecting the 
hot and cold regions. From Figure |8]^a) we see that case (1) shows much slower heat transfer 
rates compared to what is seen in case (2). This results because reconnected fleld lines in 
case (2) are linked to the globally imposed background fleld that in turn linking the hot and 
cold reservoirs. In case (1) reconnection only leads to larger fleld loops but cannot provide 
pathways between the reservoirs. The effect of different scale lengths on the evolving fleld 
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loops is shown in Figure [8[b) in which we plot the resuh from three simulations wavelengths 
for the tangled field component (tangled field "loops"). Note that A is defined in Eq (12) 
and Eq.(13). We use a sequence of values for wavelength: 2A, A and A/2. Figure [sj^b) 
clearly shows that smaller field loop A leads to the largest average heat flux, since smaller 
scale loops will reconnect before large loops for a given magnetic resistivity. This result 
demonstrates the link between the number of reconnection sites of the field and heat flux. 

We next analyze the temperature equilibration in detail. The averaged temperature 
difference across the interface is plotted in Figure [Tf^b). It shows the difference between 
the averaged temperature at the hot side and the cold side. One significant feature in 
Figure ^h) is that the temperature difference decreases to a steady value Tend in all cases. 
This resembles percolation across a membrane which allows a density jump to happen 
when filtering two fluids. Figure shows the distance required for the temperature to 
drop 80 percent at the interface. This distance characterizes the length of the interaction 
region. Except for the vertical field case where no heat transfer is allowed, the interface is 
expanding at different rates for different R values. The expansion for all the cases of nonzero 
R approaches a steady value which is also a characteristic feature of the temperature 
equilibration evolution. 

We now analyze the modification of magnetic field configuration during the evolution. 
Throughout our simulations, the local magnetic field is initially a set of complete loops 
surrounding the interaction region. Once the energy transfer begins, the interaction region 
tends to expand as discussed previously. This expansion stretches the field lines on the x 
direction and distorts these circular loops, eventually inducing magnetic reconnection which 
oppens up channels connecting the hot and cold regions. From the current Jb = |V x B|, 
we can get information on how tangled the field is. Figure ^d) shows the evolution of the 
mean value of the strength of V x B in the interaction region. We observe that in the 
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vertical and straight field case, |V x B| remains constant, but decreases to a fixed value for 
R > 2 cases. This means the field in high R cases is straightened by the stretching of the 
interaction region as seen in Figure |7[c). For the R < 1 cases, we see that |V x B| increases. 
This rise is due to magnetic energy brought in via the cold mass fiow and the creation of 
fine field structures that amplify Jb faster than dissipation caused by interface expansion. 

EDITOR: PLACE FIGURE [9] HERE. 

The local field distortion can be clearly demonstrated by studying the energy evolution 
of magnetic energy stored in different field components. In Figure |9](a), we plot the 
evolution of mean magnetic energy stored in the vertical field By/2, compared with B'^/2. 
We note that the latter includes only the fiuctuating contribution to the energy in the x 
field- that is, the contribution to the horizontal field that does not come from the global 
mean x component. 

From Figure |9]^a), we observe that the By energy decreases while the energy either 
increases or remains the same for all cases. The magnetic energy evolution can thus be 
viewed as a conversion of vertical field to horizontal field. This conversion need not conserve 
the total magnetic energy of the local tangled field because of magnetic reconnection and 
because material advecting magnetic field can fiow in and out of the interaction region. By 
comparison, in the R > 1 cases, the thermal energy and local magnetic energy can both 
decrease and add to the kinetic energy of the material surrounding the interface, because of 
the fast thermal diffusion enabled by the strong global field. 

The distortion of the local field loops can also be demonstrated by plotting the mean 
eccentricity of the field loops. In Figure [9|^b), we plot the mean eccentricity evolution. For 
all cases, the mean eccentricity is zero initially because of the circular shape of the field 
loops. Later in the evolution, large R cases tend to evolve into a state of large eccentricity 
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in the steady state. This is caused by a rapid expansion of the interface induced by the 
strong global field. In short, large R induces more distorted local field loops and less tangled 
total field due to fast interface expansion, while small R values results in less eccentric local 
field loops but with more tangled total field and strong magnetic reconnection. 



EDITOR: PLACE FIGURE M HERE. 



To compute the estimated heat transfer rate in the simulation, we calculate the 
averaged slope of the curve plotted in Figure [7||^b), and compare it to the analytic model 
in Section 4. Although the equilibration rate represented by the slope of the curves in 
Figure is changing throughout the evolution, an early phase of the evolution can be 
chosen when the field configuration has not been modified significantly for which we can 
then comptute the averaged heat transfer rate. By normalizing the resulting heat transfer 
rate to the isotropic value, we can determine the heat transfer efficiency for different 



magnetic structures. From Figure 10, we can see that the analytic prediction and the 
simulation results agree quite well except for the situation when R is below 0.2. The 
simulation result does not converge to point (0,0) but ends at an intercept on the y axis. 
This intercept indicates that even if there are initially negligibly few channels for energy 
transfer, the magnetic reconnection can open up channels and allow heat transfer. Eq.(19) 
is valid for predicting the cooling rate of the hot material throughout the early phase of the 
heat equilibration process. It also provides insight on the strength of the local field in the 
vicinity of the interface once we know the cooling rate and global magnetic field strength. 

To summarize our results we find first that the average heat flux at the end of our 
simulations is lower than at the beginning for all R values. Thus we see an approach to 
thermal equilibrum. In some cases we also see that the heat shows an initially increasing 
phase denoting a period of active magnetic reconnection. 
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In the simulations we see the average temperature difference decreases to a constant 
value Tend which is related to R. We also see the width of the initial interface expand to a 
fixed value during the simulation. 

Analysis of the simulation behavior shows that Jb is an accurate measure of structural 
change in the magnetic field. Current decreases to a constant value for large R cases and 
increases to a constant value for small R values. 

Finally we have shown that Eq.(19) can be used to estimate the energy transfer rate 
for an initially complicated field structure by considering the relative strength of the local 
field and the global field. For those cases for which R approaches 0, Eq.(19) becomes 
invalid since the energy transfer in is mainly induced by a feedback from the magnetic field 
reconnection. By comparing cases with different field loop length scales, we demonstrate 
that the smaller the field loop length scale, the faster the reconnection rate. 

7. Astrophysical Applications 

EDITOR: PLACE TABLE □ HERE. 

The issue of magnetized conduction fronts and their mediation of temperature 
distributions occurs in many astrophysical contexts. One long-standing problem that 
may involve anisotropic heat conduction are hot bubble temperatures in Wind Blown 
Bubbles (WBB). WBB's occur in a number of setting including the Planetary Nebula (PN), 
Luminious Blue Variables (LBVs) and environments of Wolf-Rayet stars. When a central 
source drives a fast wind (K^md ~ 500A;m/s) temperatures in the shocked wind material are 
expected to be of order 10^ K, which is greater than 2 kev. The temperatures observed in 
many WBB hot bubbles via from X ray emission are, however in the range of 0.5 kev to 
1 kev range. NGC 6888 is a particularly well known and well studied example for a WR star 
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(Zhekov & Park (2010)). For planetary nebulae, Chandra X-ray observations have found a 



number of WBB hot bubbles with temperatures lower than expected based on fast wind 



speeds (Montez et al. (2005), Kastner et al. (2008)). The role of wind properties and heat 
conduction in reducing hot bubble temperatures has been discussed by a number of authors 



(Steffen et al. (2008), Akashi et al. (2007), Stute fc Sahai (2007)). The role of magnetic 



fields and heat conduction was discussed in Soker (1994). 



While our simulations herein were meant to be idealized experiments aimed at 
identifying basic principles of anisotropic heat conduction fronts, we can apply physical 
scales to the simulations in order to make contact with WBB evolution. Tab.l shows the 
results of such scaling. Upon doing so, we infer that: (1) given field strengths expected for 
WBB's, heat conduction is likely to be strong enough to influence on the temperature of 
the expanding hot bubble and the cold shell bounding it. We also note that magnetic fields 



in WBB (for PN field strengths see Wouter et al. (2006)) are usually in the milli-Gauss 



range, and are relatively much stronger than the field strength that can be scaled to 
our simulations. Thus the magnetic field in realistic WBBs is highly likely to result in 
anisotropicity and regulate the behavior of heat conduction. Since the heat transfer does 
not directly depend on the magnetic /3, we can thus apply our analysis to the WBB interface 
if we approximate the interface to be planar and stationary, which is reasonable as the 
radius of curvature of WBBs are much greater than the interface scale of relevance. We 
must also assume that the global magnetic field is primarily radial. 

The computational parameters used in our simulations and the real physics parameters 
typical in a WBB are listed in the first two columns of Tab.l. We choose the domain length 
to be 0.025 pc, which is about 1 percent of the radius of the actual WBB. Tab.l shows 
that by choosing the proper scaling, our simulation fits well with the data observed in a 
typical WBB. Therefore, the conclusions we draw by analysing the simulation results and 
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the analytical expressions, especially Eq.(19), can be helpful in analyzing WBB evolution. 

8. Conclusion 

We have investigated the problem of heat transfer in regions of initially arbitrarily 
tangled magnetic fields in laminar high f3 MHD fiows using simulation results of 
ASTROBEAR code with anisotropic heat conduction. Three conclusions stand out: 

(1) Hot and cold regions initially separated by a tangled field region with locally 
confined field loops may still evolve to incur heat transfer. The local redistribution of fiuid 
elements bend the field lines and lead to magnetic reconnection that can eventually connect 
the hot and cold regions on the two sides. (2) The temperature gradient through such a 
penetrated tangeld field region tends to reach a steady state that depends on the energy 
difference between the hot and cold reservoirs on the two ends. (3) Eq.(19), a measure of 
the initial field tangle, is a good predictor of the ultimate heat transfer efficiencies across 
the interface for a wide range of R. 

A basic limitation of our simulations is that they are 2-D. A 3-D version of this study 
would be of interest as the field would then have finite scales in the third dimension possibly 
allowing channels for heat transfer excluded in 2-D. We have also not considered the effects 
of cooling in our simulations. 

Future directions of analysis could include a multi-mode study, which investigates 
the effect of the spatial spectrum of the magnetic field distribution on the heat transfer 
efficiency When there are multiple modes or a spectrum is continuous, it would be useful to 
predict how the efficiency would depend on the spectrum. In this context, a more detailed 
comparison of heat transfer in initially laminar versus initially turbulent systems would be 
of interest. 
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Appendix: Code Test 

The MHD solver and the linear thermal diffusion solver are verified by well-known tests 
such as the field loop convection problem and the Guassian diffusion problem separately. 
As a comprehensive test that involves both MHD and thermal diffusion, we use the 
magneto-thermal instability (MTI) problem to test the accuracy of the ASTROBEAR 



code with anisotropic heat conduction (Parrish & Stone (2005), Cunningham et al. (2009)). 



The problem involves setting up a 2-D temperature profile with uniform gravity pointing 
on the y direction. The domain is square with length of 0.1 in compuational units. The 
temperature and density profiles are: 

T = T,{l-y/y,) (1) 

p = Po (1 - y/yo? (2) 
EDITOR: PLACE FIGURE [iT] HERE. 



EDITOR: PLACE FIGURE [12] HERE. 



EDITOR: PLACE FIGURE M HERE. 



with i/Q = 3. The pressure profile is set up so that a hydrostatic balance may be 
achieved with uniform gravity with gravitational acceleration g = 1 in computational units. 
We also set Tq = 1 and po = 1 in computational units. There is a uniform magnetic field on 
the X direction with field strength i?o = 1-0 x 10^'^ in computational units. The anisotropic 
heat conductivity is set to be k = 1 x 10~^ in computational units. We use the pressure 
equilibrium condition for the top and bottom boundaries, that is, the pressure in the ghost 
cells are set so that its gradient balances the gravitational force. On the x direction, we use 
the periodic boundary condition. 

Initially, the domain is in pressure equilibrium. We then seed a small velocity 
perturbation: 



Vper = vq sin(n7rx/A) (3) 

with f = 1 X 10^^ and A = 0.5. This perturbation will cause the fiuid elements to have 
a tiny oscillation on y axis as well as the field lines. Once the field lines are slightly bent, 
they open up channels for heat to transfer on the y direction thus allowing the heat on the 
lower half of the domain to fiow to the upper half. It can be shown that this process has a 
positive feedback so that once the heat exchange happens, more channels will be openned 
up for heat conduction. Therefore this process forms an instability whose growth rate can 
be verified according to the linear instability growth theory. We use to denote the sound 



crossing time for the initial state. Figure 11 shows the time evolution of the field lines at 
various stages in our MTI simulation. 

We study the MTI growth rate by considering the acceleration of the fiuid elements. 
The mean speed on the y direction for the fiuid should follow the exponential growth: 
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where Vper is the strength of the initial velocity perturbation applied, 7 denotes the 
growth rate in the linear regime. We obtain the growth rate 7 by plotting In Vy against 
the evolution time and then measuring the local slope through a certain time span. The 
In Vy vs t curve is plotted in Figure [l2|^a), which shows a nice linear relation. We plot the 
growth rate against evolution time. It should be stable around the theoretical value 0.4 



initially and then decrease sharply due to the nonlinear effect. Figure 12 b). shows that the 
simulation meets our expectation fairly well. 

We also look at the energy evolution in the linear regime. The mean kinetic energy 
should first stay stable and then enter into an exponential growing phase until it hits a 
cap at around t = 200 which denotes the starting of the nonlinear phase. The evolution of 
magnetic energy should follow similar pattern as to the kinetic energy evolution, but lagged 



behind. In Figure 13, we plot the time evolution of the mean kinetic and magnetic energy 



evolutions. The results confirms the physical intuition quite well. 
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(a) 




Fig. 1. — The initial and steady state field configuration, (a): the initial field forms complete 
loops that only allows heat transfer within the interaction region, (b): the steady state field 
reconnects itself so that it allows heat transfer between regions deeply into the hot and cold 
areas. 
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Fig. 2. — Evolution of temperature distribution with R = 0.0. The cuts are at (a): t = 0.0, 
the initial state, (b): t = 0.4, (c): t = 0.8, (d): t = 1.2, the steady state. 




Fig. 3. — Evolution of density distribution with R = 0.0. The cuts are at (a): t = 0.0, the 
initial state, (b): t = 0.4, (c): t = 0.8, (d): t = 1.2, the steady state. 
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Temperature 





Fig. 4. — Evolution of temperature distribution with R = 0.4. The cuts are at (a): t = 0.0, 
the initial state, (b): t = 0.4, (c): t = 0.8, (d): t = 1.2, the steady state. 
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Fig. 5. — Evolution of temperature distribution with R = 1.0. The cuts are at (a): t = 0.0, 
the initial state, (b): t = 0.4, (c): t = 0.8, (d): t = 1.2, the steady state. 
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Fig. 6. — Evolution of mean cut temperature averaged on y direction with different R values 
labeled by different colors. The cuts are at (a): t = 0.0, the initial state, (b): t = 0.4, (c): 
t = 0.8, (d): t = 1.2, the steady state. 
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Fig. 7. — (a) top left: time evolution of mean heat flux at the interface, (b) top right: time 
evolution of average temperature difference between the hot and cold regions, (c) bottom 
left: time evolution of interface width, (d) bottom right: time evolution of the mean value 
of IV X B|. 




Fig. 8. — (a) Comparison of averaged heat flux for situation with field loops filling up the 
entire domain and situation with field loops only fill the interaction region, (b) Comparison 
of averaged heat flux for situations with different tangeld field length scale. 
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Fig. 9. — (a) Comparison on evolution of local field energy in terms of B^. and By. Circles 
corresponds to the B'^/2 curve, stars corresponds to the By/2 curve. The different colors 
denote various R values, (b) Eccentricity of the ellipses constructed by assigning the mean 
values of local \Bx\ and \By\ to the major and minor axes, respectively. The set of curves 
show different evolution patterns for different R values. 




Fig. 10. — heat transfer rate observed in the simulation compared with the analytic model. 
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Fig. 11. — Field line evolution of magneto-thermal instability, (a): initial state, (b) 
75r«. {c)t = 150r,. (d):t = 250r,. 



- 39 - 



InV vs Time 




Growth Rate werTime 




180 200 



(b) 



Fig. 12. — (a): In Vy against evolution time in r^. (b): calculated growth rate against evolu- 
tion time in computational units. Initially the growth rate is stable around the theoretical 
value 0.4 and then decreases sharply after t = 200, which indicates the evolution has entered 
the nonlinear regime. 




Fig. 13. — (a): evolution of mean kinetic energy, (b): evolution of mean magnetic energy. 
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Tabic 1: Scaling of Simulation Parameters 



Variables 


Computional Units 


WBB 


Number Density 


1 


1 cm-3 


Temperature 


100 


1 kev 


Domain Length 


0.1 


0.025 pc 


Local Field Strength 


10-3 


2-^ Gauss 


Global Field Strength 


10-^ 


Gauss 


Evolution Time 


1.2 


12,000 yrs 


Heat Conductivity 


10-2 


2 X 10-^^ cm s g-^ K-^-^ 



